#Democratic Reputations in Crisis and War: replication data
#November 20, 2021
#Josh Kertzer

######### Block A Replication III 112021.R: this R file contains the code necessary to replicate the analysis in the supplementary appendix - run Block A Replication I 112021.R and Block A Replication II 112021.R first!
# All of the following analyses were carried out using R version 4.1.2 GUI 1.77 (Big Sur ARM build 8007) on an M1 Pro Macbook Pro running MacOS Monterey. 

##### Appendix C

## Table 1: Regression model: Knesset sample

mod.1a <- lm(DV_A1 ~ A, data=knes)
mod.1b <- lm(DV_A1 ~ A + FAexp1 + milAssert1 + intTrust1 + activeCombat + ideo1 + hawk, data=knes)
mod.2a <- lm(DV_A2 ~ A, data=knes)
mod.2b <- lm(DV_A2 ~ A + FAexp1 + milAssert1 + intTrust1 + activeCombat + ideo1 + hawk, data=knes)

stargazer(mod.1a, mod.1b, mod.2a, mod.2b, omit.stat=c("LL", "ser", "f"), style="apsr", digits=3, label="tab:knessetresults", covariate.labels=c("Democracy", "Foreign affairs experience", "Military assertiveness", "International trust", "Combat experience", "L-R ideology", "Hawk (Arab/Israeli conflict)", "Constant"), dep.var.labels.include=FALSE, column.labels=c("Resolve", "Military Effectiveness"), column.separate=c(2,2))

### Table 2: Israeli Public Samples


#Male:
tab2_1 <- round(c(table(isr1$Male)[2]/sum(table(isr1$Male)), table(isr2$Male)[2]/sum(table(isr2$Male))),digits=3) #53.2%, 52.4%

#Education
tab2_2 <- round(cbind(table(isr1$educ)/sum(table(isr1$educ)), table(isr2$educ)/sum(table(isr2$educ))),digits=3) #2.1%, 1.7%; 33.2%, 25.8%; 22.5%, 22.3%; 27.0%, 31.1%; 13.6%, 17.3%; 1.6%, 1.9%

#Military experience
tab2_3 <- round(cbind(table(isr1$army)/sum(table(isr1$army)), table(isr2$army)/sum(table(isr2$army))),digits=3) #21.6%, 19.7; 50.3%, 49.7%; 28.1%, 30.7%

#Location
tab2_4 <- round(cbind(table(isr1$residence)/sum(table(isr1$residence)), table(isr2$residence)/sum(table(isr2$residence))),digits=3) #9.6%, 9.9%; 12.4%, 14.6%; 18.0 %, 17.9%; 15.1%, 13.6%; 11.2%, 13.0%; 12.1%, 8.7%; 10.2%, 10.9%; 8.3%, 9.1%; 3.0%, 2.3%

#Religion
tab2_5 <- round(cbind(table(isr1$religion)/sum(table(isr1$religion)), table(isr2$religion)/sum(table(isr2$religion))),digits=3) #60.4%, 53.5%; 19.4%, 29.9%; 15.1%, 14.0%; 5.1%, 2.6%

#Birth country
tab2_6 <- round(cbind(table(isr1$brthCtry)/sum(table(isr1$brthCtry)), table(isr2$brthCtry)/sum(table(isr2$brthCtry))),digits=3) #81.3%, 80.6%; 9.6%, 10.5%; 9.1%, 8.9%

#Means:
#Age
tab2_7 <- c(round(mean(isr1$age,na.rm=TRUE),digits=2),round(mean(isr2$age,na.rm=TRUE),digits=2)) #41.8, 41.8
#Military assertiveness
tab2_8 <- c(round(mean(isr1$milAssert1,na.rm=TRUE),digits=3), round(mean(isr2$milAssert1,na.rm=TRUE),digits=3)) #0.56, 0.58
#Right wing ideology
tab2_9 <- c(round(mean(isr1$ideo1,na.rm=TRUE),digits=3), round(mean(isr2$ideo1,na.rm=TRUE),digits=3)) #0.59, 0.61
#Hardline (Arab-Israeli)
tab2_10 <- c(round(mean(isr1$hawk,na.rm=TRUE),digits=3), round(mean(isr2$hawk,na.rm=TRUE),digits=3)) #0.62, 0.63
#International trust
tab2_11 <- c(round(mean(isr1$intTrust1,na.rm=TRUE),digits=3),NA) #0.32

tab2 <- rbind(tab2_1, tab2_2, tab2_3, tab2_4, tab2_5, tab2_6, tab2_7, tab2_8, tab2_9, tab2_10, tab2_11)
colnames(tab2) <- c("Israel I", "Israel II")
rownames(tab2) <- c("Male", "No High School degree", "High School degree", "Some college", "College degree", "Masters degree", "Doctoral degree", "Did not serve", "Served, no active combat", "Combat experience", "Jerusalem", "Tel Aviv", "Central Zone", "Haifa", "Northern Region", "Southern Region", "Lowland", "Sharon Area", "Yehuda and Shomron", "Secular", "Traditional", "Religious", "Orthodox", "Israel", "Former USSR", "Other", "Age", "Military Assertiveness", "Right Wing Ideology", "Hawkishness (Arab-Israeli conflict)", "International Trust")
stargazer(tab2, digits=2)

#SDs:
#Age
c(round(sd(isr1$age,na.rm=TRUE),digits=2), round(sd(isr2$age,na.rm=TRUE),digits=1)) #14.6,  #14.7
#Military assertiveness
c(round(sd(isr1$milAssert1,na.rm=TRUE),digits=2), round(sd(isr2$milAssert1,na.rm=TRUE),digits=2)) #0.19 #0.20
#Right wing ideology
c(round(sd(isr1$ideo1,na.rm=TRUE),digits=2), round(sd(isr2$ideo1,na.rm=TRUE),digits=2)) #0.25, 0.23
#Hardline (Arab-Israeli)
c(round(sd(isr1$hawk,na.rm=TRUE),digits=2), round(sd(isr2$hawk,na.rm=TRUE),digits=2)) #0.25  #0.25, 
#International trust
round(sd(isr1$intTrust1,na.rm=TRUE),digits=2) #0.27

#### Table 3: Elite-Mass Comparison

tab3_1 <- round(c(mean(knes$milAssert1[which(knes$currentMK==1)],na.rm=TRUE), mean(knes$ideo1[which(knes$currentMK==1)], na.rm=TRUE), mean(knes$hawk[which(knes$currentMK==1)], na.rm=TRUE), mean(knes$intTrust1[which(knes$currentMK==1)], na.rm=TRUE)),digits=2)
tab3_2 <- round(c(mean(knes$milAssert1[which(knes$currentMK==0)],na.rm=TRUE), mean(knes$ideo1[which(knes$currentMK==0)], na.rm=TRUE), mean(knes$hawk[which(knes$currentMK==0)], na.rm=TRUE), mean(knes$intTrust1[which(knes$currentMK==0)], na.rm=TRUE)),digits=2)
tab3_3 <- round(c(mean(knes$milAssert1,na.rm=TRUE), mean(knes$ideo1, na.rm=TRUE), mean(knes$hawk, na.rm=TRUE), mean(knes$intTrust1, na.rm=TRUE)),digits=2)
tab3_4 <- round(c(mean(isr1$milAssert1,na.rm=TRUE), mean(isr1$ideo1, na.rm=TRUE), mean(isr1$hawk, na.rm=TRUE), mean(isr1$intTrust1, na.rm=TRUE)),digits=2)
tab3_5 <- round(c(mean(isr2$milAssert1,na.rm=TRUE), mean(isr2$ideo1, na.rm=TRUE), mean(isr2$hawk, na.rm=TRUE), mean(isr2$intTrust1, na.rm=TRUE)),digits=2)
tab3_6 <- tab3_3-tab3_4
tab3_7 <- tab3_3-tab3_5

tab3 <- cbind(tab3_1, tab3_2, tab3_3, tab3_4, tab3_5, tab3_6, tab3_7)
rownames(tab3) <- c("Military Assertiveness", "Right Wing Ideology", "Hawkishness (Arab/Israeli)", "International Trust")

stargazer(tab3)

#Use the rank-sum tests below to calculate significant differences
#Elites vs public 1
knes.min <- cbind(with(knes, cbind(milAssert1,ideo1,hawk,intTrust1)),elite=1)
pub.min <- cbind(with(isr1, cbind(milAssert1,ideo1,hawk,intTrust1)),elite=0)
knes.pub1 <- rbind(knes.min, pub.min)

wilcox.test(milAssert1 ~ elite, data=knes.pub1) # p < 0.11
wilcox.test(ideo1 ~ elite, data=knes.pub1) # p < 0.000
wilcox.test(hawk ~ elite, data=knes.pub1) # p < 0.000
wilcox.test(intTrust1 ~ elite, data=knes.pub1) # p < 0.002

#Elites vs public 2
knes2.min <- cbind(with(knes, cbind(milAssert1,ideo1,hawk)),elite=1)
pub2.min <- cbind(with(isr2, cbind(milAssert1,ideo1,hawk)),elite=0)
knes.pub2 <- rbind(knes2.min, pub2.min)

wilcox.test(milAssert1 ~ elite, data=knes.pub2) # p < 0.02
wilcox.test(ideo1 ~ elite, data=knes.pub2) # p < 0.000
wilcox.test(hawk ~ elite, data=knes.pub2) # p < 0.000

#### Table 4: Regression models: Israeli public sample

mod.1a <- lm(DV_A1 ~ A, data=isr1)
mod.1b <- lm(DV_A1 ~ A + Male + age + educ1 + know1 + combat + milAssert1 + intTrust1 + ideo1 + hawk + isrlBrn, data=isr1)
mod.1c <- lm(DV_A1 ~ A, data=isr2)
mod.1d <- lm(DV_A1 ~ A + Male + age + educ1 + know1 + combat + milAssert1 + ideo1 + hawk + isrlBrn, data=isr2)
mod.2a <- lm(DV_A2 ~ A, data=isr1)
mod.2b <- lm(DV_A2 ~ A + Male + age + educ1 + know1 + combat + milAssert1 + intTrust1 + ideo1 + hawk + isrlBrn, data=isr1)
mod.2c <- lm(DV_A2 ~ A, data=isr2)
mod.2d <- lm(DV_A2 ~ A + Male + age + educ1 + know1 + combat + milAssert1 + ideo1 + hawk + isrlBrn, data=isr2)
stargazer(mod.1a, mod.1b, mod.1c, mod.1d, mod.2a, mod.2b, mod.2c, mod.2d, omit.stat=c("LL", "ser", "f"), style="apsr", digits=3, label="tab:publicresults", covariate.labels=c("Democracy", "Male", "Age", "Education", "Political knowledge", "Combat experience", "Military assertiveness", "International trust", "Ideology", "Hawk (Arab-Israeli)", "Born in Israel", "Constant"))

#### Appendix C.4: Heterogeneous effects in the Israeli sample

#Footnote 3: Region in which study took place

isr2$Asia <- as.numeric(isr2$gRegion_1)
isr2$America <- as.numeric(isr2$gRegion_2)
isr2$Europe <- as.numeric(isr2$gRegion_3)
isr2$Africa <- as.numeric(isr2$gRegion_4)
isr2$MENA <- as.numeric(isr2$gRegion_5)

#Yilmaz's correlation coefficient
apcorr.nosampling <- function( truth, estimate )
{
  # we're going to say that the scores at rank i are for an
  # item with an ID of i.  Thus "ID" means the index into
  # the truth and estimate vectors.
  n <- length(truth) 
  if ( length(estimate) != n )
     stop( "must be same length" )
  truth.order <- order( truth, decreasing=TRUE ) 
  estimate.order <- order( estimate, decreasing=TRUE ) 
  innerSum <- 0
  for ( i in 2:n )
  {
	currDocID <- estimate.order[i] 
	estimate.rankedHigherIDs <- estimate.order[1:(i-1)] 
	# where is the current doc in the truth order?
	currDoc.truth.order.index <- which( truth.order == currDocID )
	truth.rankedHigherIDs <- vector()
	if ( currDoc.truth.order.index != 1 ) # top ranked doc, beware
	{
	    truth.rankedHigherIDs <- truth.order[1:(currDoc.truth.order.index-1)]
	}
	C_i <- length( intersect(estimate.rankedHigherIDs, truth.rankedHigherIDs) )
	innerSum <- innerSum + (C_i / (i-1))
  }
  result = 2 / (n-1) * innerSum - 1   
  return( result )
}

#Yilmaz's correlation coefficient
rank.dat <- 6+cbind(isr2$Asia, isr2$America, isr2$Europe, isr2$Africa, isr2$MENA)*-1 #Reverse it: higher values = higher rankings
apply(rank.dat,2,mean) #Let's set the "true" ranking to the rank of the average rank for each region: 3,5,4,2,1 (i.e., Middle East is the most popular, America the least) --> 3,1,2,4,5
yilmaz <- rep(NA, nrow(rank.dat))
kendall <- rep(NA, nrow(rank.dat))
for (i in 1:nrow(rank.dat)){
	yilmaz[i] <- apcorr.nosampling(truth=c(3,1,2,4,5),estimate=rank.dat[i,])
	kendall[i] <- cor.test(x=c(3,1,2,4,5), y=rank.dat[i,], method="kendall")$estimate
}

isr2$countryRank <- yilmaz
isr2$countryRank2 <- kendall

mod.1b <- lm(DV_A1 ~ A*countryRank+age + Male + educ1 + ideo1, data=isr2)
summary(mod.1b) #p < 0.85

mod.2b <- lm(DV_A2 ~ A*countryRank+age + Male + educ1 + ideo1, data=isr2)
summary(mod.2b) #p < 0.70

mod.1b <- lm(DV_A1 ~ A*countryRank2+age + Male + educ1 + ideo1, data=isr2)
summary(mod.1b) #p < 0.59

mod.2b <- lm(DV_A2 ~ A*countryRank2+age + Male + educ1 + ideo1, data=isr2)
summary(mod.2b) #p < 0.49

###### Appendix C.5: Downsampling Israeli public results

set.seed(43215)
downsample1 <- rep(NA, B)
for (i in 1:B){
  j <- sample(1:nrow(isr1),88, replace=TRUE)
  temp <- isr1[j,]
  downsample1[i] <- mean(temp$DV_A1[temp$A==1], na.rm=TRUE) - mean(temp$DV_A1[temp$A==0], na.rm=TRUE) 
}

set.seed(43215)
downsample2 <- rep(NA,B)
for (i in 1:B){
	 j <- sample(1:nrow(isr2), 88, replace=TRUE)
	 temp <- isr2[j,]
 	downsample2[i] <- mean(temp$DV_A1[temp$A==1], na.rm=TRUE) - mean(temp$DV_A1[temp$A==0], na.rm=TRUE) 
}

#Compare with boot.k1 from Block A Replication II 112021.R

x_2 <- data.frame(y=c(downsample2, boot.k1, downsample1), z=rep(c("Isr Public II (Downsampled)","Knesset", "Isr Public I (Downsampled)"), each=B), dv=rep(c("In Crises..."),each=3*B))

#### Figure 2: Downsampled crisis results

dev.new(height=3,width=5)
dv3 <- ggplot(x_2, aes(x=y)) + geom_density(aes(fill=z), alpha=0.65) + labs(x="Effect of democracy", y="Density") + scale_fill_brewer(guide=guide_legend(title="Sample"), palette="Set1") + facet_wrap(~dv ,ncol=2, scales="fixed") + theme_bw() + scale_y_continuous(breaks=NULL, limits=c(0,0.5)) + geom_vline(xintercept=0, linetype=3, alpha=0.5, show.legend=FALSE) 
dv3

######## Appendix C.6: Comparison of effect size magnitude with large-N work

### Table 5: The probability of crisis reciprocation: Bilateral MIDs
#Logits
mod.3l <- glm(recip ~ democ_a + democ_b + demdem+  majmaj+ majmin+ minmaj+ capshare_a+ contig+ swt_dyad+ tauleader_a+ tauleader_b+ territory+ government+ policy+ other + dummy1914 + dummy1915 + dummy1916 + dummy1917 + dummy1918 + dummy1939 + dummy1940 + dummy1941 + dummy1942 + dummy1943 + dummy1944 + dummy1945, data=subset(downes, downes$year<=1984), family=binomial(logit))
mod.4l <- glm(recip ~ democ_a + democ_b + demdem+  majmaj+ majmin+ minmaj+ capshare_a+ contig+ swt_dyad+ tauleader_a+ tauleader_b+ territory+ government+ policy+ other, data=subset(downes, downes$year<=1984 & downes$worldwar==0),family=binomial(logit))

#LPMs
mod.3 <- lm(recip ~ democ_a + democ_b + demdem+  majmaj+ majmin+ minmaj+ capshare_a+ contig+ swt_dyad+ tauleader_a+ tauleader_b+ territory+ government+ policy+ other + dummy1939 + dummy1940 + dummy1941 + dummy1942 + dummy1943 + dummy1944 + dummy1945 + dummy1914 + dummy1915 + dummy1916 + dummy1917 + dummy1918, data=subset(downes, downes$year<=1984))
mod.4 <- lm(recip ~ democ_a + democ_b + demdem+  majmaj+ majmin+ minmaj+ capshare_a+ contig+ swt_dyad+ tauleader_a+ tauleader_b+ territory+ government+ policy+ other, data=subset(downes, downes$year<=1984 & downes$worldwar==0))


stargazer(mod.3l,mod.3, mod.4l, mod.4, type="latex", omit.stat=c("LL", "ser", "f"), style="apsr", digits=3, label="tab:downes", covariate.labels=c("Democratic Initiator", "Democratic Target", "Both Democratic", "Major Power Initiator-Major Power Target", "Major Power Initiator-Minor Power Target", "Minor Power Initiator-Major Power Target", "Initiator's share of capabilities", "Contiguous", "Alliance portfolio similarity", "Status quo evaluation of initiator", "Status quo evaluation of target", "Territory", "Government", "Policy", "Other"), omit=c("dummy1914", "dummy1915","dummy1916", "dummy1917", "dummy1918", "dummy1939", "dummy1940", "dummy1941", "dummy1942", "dummy1943", "dummy1944", "dummy1945"))


#mod.3: -6.6% (when target is dictatorship)
#mod.4: -7.9% (when target is dictatorship)

### Table 6: Predicted probabilities of crisis reciprocation

#Maj-maj dyad
x.1.m1 <- c(1,1,0,0,1,0,0,mean(downes$capshare_a[which(downes$majmaj==1)], na.rm=TRUE),1,mean(downes$swt_dyad[which(downes$majmaj==1)], na.rm=TRUE), mean(downes$tauleader_a[which(downes$majmaj==1)],na.rm=TRUE), mean(downes$tauleader_b[which(downes$majmaj==1)],na.rm=TRUE), 0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0)  #Maj-maj dyad,nondemocratic target, contiguous, policy revision
x.0.m1 <- x.1.m1
x.0.m1[2] <- 0

#Maj-min dyad
x.1.m2 <- c(1,1,0,0,0,1,0,mean(downes$capshare_a[which(downes$majmin==1)], na.rm=TRUE),1,mean(downes$swt_dyad[which(downes$majmin==1)], na.rm=TRUE), mean(downes$tauleader_a[which(downes$majmin==1)],na.rm=TRUE), mean(downes$tauleader_b[which(downes$majmin==1)],na.rm=TRUE), 0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0)  #Maj-min dyad,nondemocratic target, contiguous, policy revision
x.0.m2 <- x.1.m2
x.0.m2[2] <- 0

#Min-maj dyad
x.1.m3 <- c(1,1,0,0,0,0,1,mean(downes$capshare_a[which(downes$minmaj==1)], na.rm=TRUE),1,mean(downes$swt_dyad[which(downes$minmaj==1)], na.rm=TRUE), mean(downes$tauleader_a[which(downes$minmaj==1)],na.rm=TRUE), mean(downes$tauleader_b[which(downes$minmaj==1)],na.rm=TRUE), 0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0)  #Min-maj dyad,nondemocratic target, contiguous, policy revision
x.0.m3 <- x.1.m3
x.0.m3[2] <- 0

#Min-min dyad
x.1.m4 <- c(1,1,0,0,0,0,0,mean(downes$capshare_a[which(downes$minmaj==0 & downes$majmin==0 & downes$majmaj==0)], na.rm=TRUE),1,mean(downes$swt_dyad[which(downes$minmaj==0 & downes$majmin==0 & downes$majmaj==0)], na.rm=TRUE), mean(downes$tauleader_a[which(downes$minmaj==0 & downes$majmin==0 & downes$majmaj==0)],na.rm=TRUE), mean(downes$tauleader_b[which(downes$minmaj==0 & downes$majmin==0 & downes$majmaj==0)],na.rm=TRUE), 0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0)  #Min-maj dyad,nondemocratic target, contiguous, policy revision
x.0.m4 <- x.1.m4
x.0.m4[2] <- 0

### World wars included

#Major-maj dyad
round(plogis(t(as.matrix(x.1.m1)) %*% as.matrix(coef(mod.3l))),digits=3) #24.6% if democratic initiator
round(plogis(t(as.matrix(x.0.m1)) %*% as.matrix(coef(mod.3l))),digits=3) #30.8% if nondemocratic initiator
round(plogis(t(as.matrix(x.1.m1)) %*% as.matrix(coef(mod.3l)))-plogis(t(as.matrix(x.0.m1)) %*% as.matrix(coef(mod.3l))),digits=3) #-6.3%

#Major-minor dyad
round(plogis(t(as.matrix(x.1.m2)) %*% as.matrix(coef(mod.3l))),digits=3) #25.4% if democratic initiator
round(plogis(t(as.matrix(x.0.m2)) %*% as.matrix(coef(mod.3l))),digits=3) #31.7% if nondemocratic initiator
round(plogis(t(as.matrix(x.1.m2)) %*% as.matrix(coef(mod.3l)))-plogis(t(as.matrix(x.0.m2)) %*% as.matrix(coef(mod.3l))),digits=3) #-6.4%

#Minor-major dyad
round(plogis(t(as.matrix(x.1.m3)) %*% as.matrix(coef(mod.3l))),digits=3) #32.5% if democratic initiator
round(plogis(t(as.matrix(x.0.m3)) %*% as.matrix(coef(mod.3l))),digits=3) #39.7% if nondemocratic initiator
round(plogis(t(as.matrix(x.1.m3)) %*% as.matrix(coef(mod.3l)))-plogis(t(as.matrix(x.0.m3)) %*% as.matrix(coef(mod.3l))),digits=3) #-7.2%

#Minor-minor dyad
round(plogis(t(as.matrix(x.1.m4)) %*% as.matrix(coef(mod.3l))),digits=3) #31.9% if democratic initiator
round(plogis(t(as.matrix(x.0.m4)) %*% as.matrix(coef(mod.3l))),digits=3) #39.1% if nondemocratic initiator
round(plogis(t(as.matrix(x.1.m4)) %*% as.matrix(coef(mod.3l)))-plogis(t(as.matrix(x.0.m4)) %*% as.matrix(coef(mod.3l))),digits=3) #-7.2%

### Now, world wars excluded

#Major-maj dyad
round(plogis(t(as.matrix(x.1.m1)[1:16]) %*% as.matrix(coef(mod.4l))),digits=3) #20.9% if democratic initiator
round(plogis(t(as.matrix(x.0.m1)[1:16]) %*% as.matrix(coef(mod.4l))),digits=3) #27.8% if nondemocratic initiator
round(plogis(t(as.matrix(x.1.m1)[1:16]) %*% as.matrix(coef(mod.4l)))-plogis(t(as.matrix(x.0.m1)[1:16]) %*% as.matrix(coef(mod.4l))),digits=3) #-6.9%

#Major-minor dyad
round(plogis(t(as.matrix(x.1.m2)[1:16]) %*% as.matrix(coef(mod.4l))),digits=3) #24.7% if democratic initiator
round(plogis(t(as.matrix(x.0.m2)[1:16]) %*% as.matrix(coef(mod.4l))),digits=3) #32.2% if nondemocratic initiator
round(plogis(t(as.matrix(x.1.m2)[1:16]) %*% as.matrix(coef(mod.4l)))-plogis(t(as.matrix(x.0.m2)[1:16]) %*% as.matrix(coef(mod.4l))),digits=3) #-7.6%

#Minor-major dyad
round(plogis(t(as.matrix(x.1.m3)[1:16]) %*% as.matrix(coef(mod.4l))),digits=3) #31.4% if democratic initiator
round(plogis(t(as.matrix(x.0.m3)[1:16]) %*% as.matrix(coef(mod.4l))),digits=3) #39.9% if nondemocratic initiator
round(plogis(t(as.matrix(x.1.m3)[1:16]) %*% as.matrix(coef(mod.4l)))-plogis(t(as.matrix(x.0.m3)[1:16]) %*% as.matrix(coef(mod.4l))),digits=3) #-8.6%

#Minor-minor dyad
round(plogis(t(as.matrix(x.1.m4)[1:16]) %*% as.matrix(coef(mod.4l))),digits=3) #30.8% if democratic initiator
round(plogis(t(as.matrix(x.0.m4)[1:16]) %*% as.matrix(coef(mod.4l))),digits=3) #39.3% if nondemocratic initiator
round(plogis(t(as.matrix(x.1.m4)[1:16]) %*% as.matrix(coef(mod.4l)))-plogis(t(as.matrix(x.0.m4)[1:16]) %*% as.matrix(coef(mod.4l))),digits=3) #-8.5%


### First, generate dichotomous variables for mosaic plots

#Knesset sample
knes$Crisis <- ifelse(knes$DV_A1 <50,"Back\nDown",ifelse(knes$DV_A1 >50,"Stand\nFirm",NA))
knes$War <- ifelse(knes$DV_A2 <50,"Lose",ifelse(knes$DV_A2 >50,"Win",NA))

#Israel I sample
isr1$Crisis <- ifelse(isr1$DV_A1 <50,"Back\nDown",ifelse(isr1$DV_A1 >50,"Stand\nFirm",NA))
isr1$War <- ifelse(isr1$DV_A2 <50,"Lose",ifelse(isr1$DV_A2 >50,"Win",NA))

#Israel II sample
isr2$Crisis <- ifelse(isr2$DV_A1 <50,"Back\nDown",ifelse(isr2$DV_A1 >50,"Stand\nFirm",NA))
isr2$War <- ifelse(isr2$DV_A2 <50,"Lose",ifelse(isr2$DV_A2 >50,"Win",NA))

#ROK sample
rok$Crisis <- ifelse(rok$DV_A1 <50,"Back\nDown",ifelse(rok$DV_A1 >50,"Stand\nFirm",NA))
rok$War <- ifelse(rok$DV_A2 <50,"Lose",ifelse(rok$DV_A2 >50,"Win",NA))

#UK sample I
uk1$Crisis <- ifelse(uk1$DV_A1 <50,"Back\nDown",ifelse(uk1$DV_A1 >50,"Stand\nFirm",NA))
uk1$War <- ifelse(uk1$DV_A2 <50,"Lose",ifelse(uk1$DV_A2 >50,"Win",NA))

#UK sample II
uk2$Crisis <- ifelse(uk2$DV_A1 <50,"Back\nDown",ifelse(uk2$DV_A1 >50,"Stand\nFirm",NA))
uk2$War <- ifelse(uk2$DV_A2 <50,"Lose",ifelse(uk2$DV_A2 >50,"Win",NA))

#USA sample
usa$Crisis <- ifelse(usa$DV_A1 <50,"Back\nDown",ifelse(usa$DV_A1 >50,"Stand\nFirm",NA))
usa$War <- ifelse(usa$DV_A2 <50,"Lose",ifelse(usa$DV_A2 >50,"Win",NA))

### Figure 3: Exploiting the two-stage structure of the experiment
## Mosaic plot panel

mosaic <- data.frame(Crisis=c(knes$Crisis, isr1$Crisis, isr2$Crisis, rok$Crisis, uk1$Crisis, uk2$Crisis, usa$Crisis), War=c(knes$War, isr1$War, isr2$War, rok$War, uk1$War, uk2$War, usa$War), A=c(knes$A, isr1$A, isr2$A, rok$A, uk1$A, uk2$A, usa$A), sample=c(rep("Israel: Knesset", nrow(knes)), rep("Israel I", nrow(isr1)), rep("Israel II", nrow(isr2)), rep("Republic of Korea", nrow(rok)), rep("United Kingdom I", nrow(uk1)), rep("United Kingdom II", nrow(uk2)), rep("USA", nrow(usa))))

#Drop observations in isr1 where regime type is missing
mosaic <- mosaic[-which(is.na(mosaic$A)),]
mosaic$A[which(mosaic$A==0)] <- "Dictatorship"
mosaic$A[which(mosaic$A==1)] <- "Democracy"

dev.new(height=3.2, width=10.2)
ggplot(data=mosaic) + geom_mosaic(aes(x=product(War, Crisis), fill=A), na.rm=TRUE, offset=0.02) + facet_grid(A ~ sample) + theme_bw() + labs(x="Crisis", y="War") + scale_x_productlist(labels=c("Back\nDown", "Stand\nFirm")) + labs(fill="Regime Type")

### Now, calculate bootstrapped product for density distribution panel

B <- 1500
set.seed(43215)

knes.mat <- matrix(NA, nrow=B, ncol=3) #Average product for democracies, dictatorship, and difference
for (i in 1:B)
{
	j <- sample(1:nrow(knes), nrow(knes), replace=TRUE)	
	temp <- knes[j,]
	knes.mat[i,1] <- 100*mean((temp$DV_A1[which(temp$A==1)]/100)*(temp$DV_A2[which(temp$A==1)]/100), na.rm=TRUE)
	knes.mat[i,2] <- 100*mean((temp$DV_A1[which(temp$A==0)]/100)*(temp$DV_A2[which(temp$A==0)]/100), na.rm=TRUE)
	knes.mat[i,3] <- 100*(mean((temp$DV_A1[which(temp$A==1)]/100)*(temp$DV_A2[which(temp$A==1)]/100), na.rm=TRUE) - mean((temp$DV_A1[which(temp$A==0)]/100)*(temp$DV_A2[which(temp$A==0)]/100), na.rm=TRUE))	
}

isr1.mat <- matrix(NA, nrow=B, ncol=3) #Average product for democracies, dictatorship, and difference
for (i in 1:B)
{
	j <- sample(1:nrow(isr1), nrow(isr1), replace=TRUE)	
	temp <- isr1[j,]
	isr1.mat[i,1] <- 100*mean((temp$DV_A1[which(temp$A==1)]/100)*(temp$DV_A2[which(temp$A==1)]/100), na.rm=TRUE)
	isr1.mat[i,2] <- 100*mean((temp$DV_A1[which(temp$A==0)]/100)*(temp$DV_A2[which(temp$A==0)]/100), na.rm=TRUE)
	isr1.mat[i,3] <- 100*(mean((temp$DV_A1[which(temp$A==1)]/100)*(temp$DV_A2[which(temp$A==1)]/100), na.rm=TRUE) - mean((temp$DV_A1[which(temp$A==0)]/100)*(temp$DV_A2[which(temp$A==0)]/100), na.rm=TRUE))	
}

isr2.mat <- matrix(NA, nrow=B, ncol=3) #Average product for democracies, dictatorship, and difference
for (i in 1:B)
{
	j <- sample(1:nrow(isr2), nrow(isr2), replace=TRUE)	
	temp <- isr2[j,]
	isr2.mat[i,1] <- 100*mean((temp$DV_A1[which(temp$A==1)]/100)*(temp$DV_A2[which(temp$A==1)]/100), na.rm=TRUE)
	isr2.mat[i,2] <- 100*mean((temp$DV_A1[which(temp$A==0)]/100)*(temp$DV_A2[which(temp$A==0)]/100), na.rm=TRUE)
	isr2.mat[i,3] <- 100*(mean((temp$DV_A1[which(temp$A==1)]/100)*(temp$DV_A2[which(temp$A==1)]/100), na.rm=TRUE) - mean((temp$DV_A1[which(temp$A==0)]/100)*(temp$DV_A2[which(temp$A==0)]/100), na.rm=TRUE))	
}

rok.mat <- matrix(NA, nrow=B, ncol=3) #Average product for democracies, dictatorship, and difference
for (i in 1:B)
{
	j <- sample(1:nrow(rok), nrow(rok), replace=TRUE)	
	temp <- rok[j,]
	rok.mat[i,1] <- 100*mean((temp$DV_A1[which(temp$A==1)]/100)*(temp$DV_A2[which(temp$A==1)]/100), na.rm=TRUE)
	rok.mat[i,2] <- 100*mean((temp$DV_A1[which(temp$A==0)]/100)*(temp$DV_A2[which(temp$A==0)]/100), na.rm=TRUE)
	rok.mat[i,3] <- 100*(mean((temp$DV_A1[which(temp$A==1)]/100)*(temp$DV_A2[which(temp$A==1)]/100), na.rm=TRUE) - mean((temp$DV_A1[which(temp$A==0)]/100)*(temp$DV_A2[which(temp$A==0)]/100), na.rm=TRUE))	
}

uk1.mat <- matrix(NA, nrow=B, ncol=3) #Average product for democracies, dictatorship, and difference
for (i in 1:B)
{
	j <- sample(1:nrow(uk1), nrow(uk1), replace=TRUE)	
	temp <- uk1[j,]
	uk1.mat[i,1] <- 100*mean((temp$DV_A1[which(temp$A==1)]/100)*(temp$DV_A2[which(temp$A==1)]/100), na.rm=TRUE)
	uk1.mat[i,2] <- 100*mean((temp$DV_A1[which(temp$A==0)]/100)*(temp$DV_A2[which(temp$A==0)]/100), na.rm=TRUE)
	uk1.mat[i,3] <- 100*(mean((temp$DV_A1[which(temp$A==1)]/100)*(temp$DV_A2[which(temp$A==1)]/100), na.rm=TRUE) - mean((temp$DV_A1[which(temp$A==0)]/100)*(temp$DV_A2[which(temp$A==0)]/100), na.rm=TRUE))	
}

uk2.mat <- matrix(NA, nrow=B, ncol=3) #Average product for democracies, dictatorship, and difference
for (i in 1:B)
{
	j <- sample(1:nrow(uk2), nrow(uk2), replace=TRUE)	
	temp <- uk2[j,]
	uk2.mat[i,1] <- 100*mean((temp$DV_A1[which(temp$A==1)]/100)*(temp$DV_A2[which(temp$A==1)]/100), na.rm=TRUE)
	uk2.mat[i,2] <- 100*mean((temp$DV_A1[which(temp$A==0)]/100)*(temp$DV_A2[which(temp$A==0)]/100), na.rm=TRUE)
	uk2.mat[i,3] <- 100*(mean((temp$DV_A1[which(temp$A==1)]/100)*(temp$DV_A2[which(temp$A==1)]/100), na.rm=TRUE) - mean((temp$DV_A1[which(temp$A==0)]/100)*(temp$DV_A2[which(temp$A==0)]/100), na.rm=TRUE))	
}

usa.mat <- matrix(NA, nrow=B, ncol=3) #Average product for democracies, dictatorship, and difference
for (i in 1:B)
{
	j <- sample(1:nrow(usa), nrow(usa), replace=TRUE)	
	temp <- usa[j,]
	usa.mat[i,1] <- 100*mean((temp$DV_A1[which(temp$A==1)]/100)*(temp$DV_A2[which(temp$A==1)]/100), na.rm=TRUE)
	usa.mat[i,2] <- 100*mean((temp$DV_A1[which(temp$A==0)]/100)*(temp$DV_A2[which(temp$A==0)]/100), na.rm=TRUE)
	usa.mat[i,3] <- 100*(mean((temp$DV_A1[which(temp$A==1)]/100)*(temp$DV_A2[which(temp$A==1)]/100), na.rm=TRUE) - mean((temp$DV_A1[which(temp$A==0)]/100)*(temp$DV_A2[which(temp$A==0)]/100), na.rm=TRUE))	
}

diff.mat <- data.frame(diff=c(knes.mat[,3], isr1.mat[,3], isr2.mat[,3], rok.mat[,3], uk1.mat[,3], uk2.mat[,3], usa.mat[,3]), sample=rep(c("Israel: Knesset", "Israel I", "Israel II", "Republic of Korea", "United Kingdom I", "United Kingdom II", "USA"), each=B))

dev.new(height=3.2, width=10.2)
ggplot(diff.mat) + geom_density(aes(x=diff, y=..density.., fill=sample), alpha=0.55) + scale_fill_brewer(guide=guide_legend(title="Sample"), palette="Set1") + labs(x="Overall effect of democracy (Pr(Stand Firm) * Pr(Win))", y="") + scale_y_continuous(breaks=NULL) + theme_bw() + geom_vline(xintercept=0, linetype=3, alpha=0.5, show.legend=FALSE)

#### Sequential g-estimation

mod.knes1 <- lm(DV_A2 ~ A, data=knes)
mod.knes1g <- sequential_g(DV_A2 ~ A | 1 | DV_A1, data=knes)

mod.isr1 <- lm(DV_A2 ~ A, data=isr1)
mod.isr1g <- sequential_g(DV_A2 ~ A | 1 | DV_A1, data=isr1)

mod.isr2 <- lm(DV_A2 ~ A, data=isr2)
mod.isr2g <- sequential_g(DV_A2 ~ A | 1 | DV_A1, data=isr2)

mod.uk1 <- lm(DV_A2 ~ A, data=uk1)
mod.uk1g <- sequential_g(DV_A2 ~ A | 1 | DV_A1, data=uk1)

mod.uk2 <- lm(DV_A2 ~ A, data=uk2)
mod.uk2g <- sequential_g(DV_A2 ~ A | 1 | DV_A1, data=uk2)

mod.rok1 <- lm(DV_A2 ~ A, data=rok)
mod.rok1g <- sequential_g(DV_A2 ~ A | 1 | DV_A1, data=rok)

mod.usa1 <- lm(DV_A2 ~ A, data=usa2)
mod.usa1g <- sequential_g(DV_A2 ~ A | 1 | DV_A1, data=usa2)


B <- c(coef(mod.isr1)[2], coef(mod.isr1g)[2], coef(mod.isr2)[2], coef(mod.isr2g)[2], coef(mod.knes1)[2], coef(mod.knes1g)[2], coef(mod.rok1)[2], coef(mod.rok1g)[2], coef(mod.uk1)[2], coef(mod.uk1g)[2], coef(mod.uk2)[2], coef(mod.uk2g)[2], coef(mod.usa1)[2], coef(mod.usa1g)[2]) 
S <- c(summary(mod.isr1)$coefficients[2,2], summary(mod.isr1g)$coefficients[2,2], summary(mod.isr2)$coefficients[2,2], summary(mod.isr2g)$coefficients[2,2], summary(mod.knes1)$coefficients[2,2], summary(mod.knes1g)$coefficients[2,2], summary(mod.rok1)$coefficients[2,2], summary(mod.rok1g)$coefficients[2,2], summary(mod.uk1)$coefficients[2,2], summary(mod.uk1g)$coefficients[2,2], summary(mod.uk2)$coefficients[2,2], summary(mod.uk2g)$coefficients[2,2], summary(mod.usa1)$coefficients[2,2], summary(mod.usa1g)$coefficients[2,2]) 

x <- rbind(B,S)
out <- as.matrix(rbind(x[,1:2], x[,3:4], x[,5:6], x[,7:8], x[,9:10], x[,11:12], x[,13:14]))
colnames(out) <- c("ATE", "ACDE")
rownames(out)[seq(1,13,by=2)] <- c("Israel (Public I)", "Israel (Public II)", "Israel (Knesset)", "ROK", "UK (Public I)", "UK (Public II)", "USA")
rownames(out)[seq(2,14,by=2)] <- ""
stargazer(out)


### Primary and secondary outcomes (not converted on probability scale)

isr2$resolve1 <- isr2$DV_A1
isr2$publicPacifist1 <- 100*rescale(isr2$publicPacifist)
isr2$elitePacifist1 <- 100*rescale(isr2$elitePacifist)
isr2$casualtySensitive1 <- 100*rescale(isr2$casualtySensitive)
isr2$finSensitive1 <- 100*rescale(isr2$finSensitive)
isr2$pubSupportDiff1 <- 50+10*(isr2$pubSupportDiff) #Within-subject scale needs to be rescaled manually
isr2$turnoverDiff1 <- 50+10*rescale(isr2$turnoverDiff) #Within-subject scale needs to be rescaled manually
isr2$leaderFate1 <- 100*rescale(isr2$leaderFate)
isr2$pubSupportDiffC1 <- isr2$pubSupportDiff1
isr2$pubSupportDiffC1[which(isr2$mech!="Closed_Crisis")] <- NA
isr2$turnoverDiffC1 <- isr2$turnoverDiff1
isr2$turnoverDiffC1[which(isr2$mech!="Closed_Crisis")] <- NA
isr2$leaderFateC1 <- isr2$leaderFate1
isr2$leaderFateC1[which(isr2$mech!="Closed_Crisis")] <- NA
isr2$followThrough1 <- 100*rescale(isr2$followThrough)

## Wars: Win, selection effects, domestic constraints (3), military conduct (3)
isr2$win1 <- isr2$DV_A2
isr2$selectWin1 <- 100*rescale(isr2$selectWin)
isr2$pubSupportDiffW1 <- isr2$pubSupportDiff1
isr2$pubSupportDiffW1[which(isr2$mech!="Closed_War")] <- NA 
isr2$turnoverDiffW1 <- isr2$turnoverDiff1
isr2$turnoverDiffW1[which(isr2$mech!="Closed_War")] <- NA
isr2$leaderFateW1 <- isr2$leaderFate
isr2$leaderFateW1[which(isr2$mech!="Closed_War")] <- NA
isr2$reliableAllies1 <- 100*rescale(isr2$reliableAllies)
isr2$training1 <- 100*rescale(isr2$training)
isr2$morale1 <- 100*rescale(isr2$morale)

c.vec <- c("resolve1","publicPacifist1","elitePacifist1","casualtySensitive1","finSensitive1", "turnoverDiffC1", "pubSupportDiffC1", "leaderFateC1", "followThrough1")
w.vec <- c("win1", "selectWin1", "turnoverDiffW1", "pubSupportDiffW1", "leaderFateW1", "reliableAllies1", "training1", "morale1")

c.vec.lab <- c("Stand firm in crises", "Dovish public", "Dovish leaders", "Casualty sensitivity", "Financial sensitivity", "Leader turnover", "Public support", "Leader fate", "Follow through")
w.vec.lab <- c("Win in war", "Pr(Victory|Initiate)", "Leader turnover",  "Public support", "Leader fate", "Allies", "Training", "Morale")


## Parallel processing (to speed it up)
cl <- makeCluster(20,"SOCK") #Modify for the number of cores available on your computer
#For bootstrapping
B <- 1500
set.seed(43215)

bootA <- function(...){
	j <- sample(1:nrow(isr2), replace=TRUE)
	temp.out <- rep(NA, length(x.vec))
	X <- isr2$A[j]
	for (k in 1:length(x.vec)){
		Y <- isr2[j,x.vec[k]]
		temp.out[k] <- mean(Y[X==1],na.rm=TRUE) - mean(Y[X==0],na.rm=TRUE)
	}
	return(temp.out)
}

clusterSetRNGStream(cl, 43215) #Set seed for parallel processing
x.vec <- c.vec
clusterExport(cl, list("isr2", "x.vec", "bootA"), envir=environment())
mat.c <- parSapply(cl,rep(1,B),bootA)

x.vec <- w.vec
clusterExport(cl, list("isr2", "x.vec", "bootA"), envir=environment())
mat.w <- parSapply(cl,rep(1,B),bootA)

mat.c2 <- t(mat.c)
colnames(mat.c2) <- c.vec.lab

mat.w2 <- t(mat.w)
colnames(mat.w2) <- w.vec.lab

df.c <- gather(as.data.frame(mat.c2))
df.c$key <- factor(df.c$key, levels=rev(c.vec.lab))

df.w <- gather(as.data.frame(mat.w2))
df.w$key <- factor(df.w$key, levels=rev(w.vec.lab))

dev.new(height=5, width=5)
ggplot(df.c, aes(y=key)) + geom_density_ridges(aes(x=value)) + theme_bw() + xlab("Effect of democracy in crises") + ylab("") + geom_vline(xintercept=0, lty=3) + xlim(-10,30)

dev.new(height=5, width=5)
ggplot(df.w, aes(y=key)) + geom_density_ridges(aes(x=value)) + theme_bw() + xlab("Effect of democracy in war") + ylab("") + geom_vline(xintercept=0, lty=3) + xlim(-10,30)

###### Appendix D.2: Manipulating Regime Type and Alliance


### Figure 5: Results with transitional democracy condition (Israeli public sample)

isr1$A2 <- isr1$A #0 = Dictatorship, 1 = Democracy
isr1$A2[which(isr1$RegimeA==3)] <- 0.5 #In between democracy and dictatorship

B <- 1500
set.seed(43215)

out.a <- matrix(NA, nrow=B, ncol=6)
for (i in 1:B){
	j <- sample(1:nrow(isr1), replace=TRUE)
	temp <- isr1[j,]
	out.a[i,1] <- mean(temp$DV_A1[which(temp$A2==0)], na.rm=TRUE)
	out.a[i,2] <- mean(temp$DV_A1[which(temp$A2==0.5)],na.rm=TRUE)
	out.a[i,3] <- mean(temp$DV_A1[which(temp$A2==1)],na.rm=TRUE)
	out.a[i,4] <- mean(temp$DV_A2[which(temp$A2==0)],na.rm=TRUE)
	out.a[i,5] <- mean(temp$DV_A2[which(temp$A2==0.5)],na.rm=TRUE)
	out.a[i,6] <- mean(temp$DV_A2[which(temp$A2==1)],na.rm=TRUE)			
}

out.a2 <- gather(as.data.frame(out.a))
out.a2$Outcome <- "Resolve in crises"
out.a2$Outcome[grep("V4|V5|V6",out.a2$key)] <- "Victory in war"
out.a2$Regime <- "Dictatorship"
out.a2$Regime[grep("V3|V6",out.a2$key)] <- "Democracy"
out.a2$Regime[grep("V2|V5",out.a2$key)] <- "Transitional Democracy"
out.a2$Regime <- as.factor(out.a2$Regime)
out.a2$Regime <- relevel(out.a2$Regime, ref="Transitional Democracy")
out.a2$Regime <- relevel(out.a2$Regime, ref="Democracy")


dev.new(height=5, width=7)
out.plot <- ggplot(out.a2, aes(y=value, fill=Regime)) + geom_density(alpha=0.5) + facet_wrap(~Outcome) + theme_bw() + ylab("Cell mean") + theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank()) + guides(fill=guide_legend(title="Regime Type"))
out.plot

#And now, in the US

set.seed(43215)
out.usa <- matrix(NA, nrow=B, ncol=6)
for (i in 1:B){
	j <- sample(1:nrow(usa), replace=TRUE)
	temp <- usa[j,]
	out.usa[i,1] <- mean(temp$DV_A1[which(temp$RegimeA=="Dictatorship")], na.rm=TRUE)
	out.usa[i,2] <- mean(temp$DV_A1[which(temp$RegimeA=="Dem. Transition")],na.rm=TRUE)
	out.usa[i,3] <- mean(temp$DV_A1[which(temp$RegimeA=="Democracy")],na.rm=TRUE)
	out.usa[i,4] <- mean(temp$DV_A2[which(temp$RegimeA=="Dictatorship")],na.rm=TRUE)
	out.usa[i,5] <- mean(temp$DV_A2[which(temp$RegimeA=="Dem. Transition")],na.rm=TRUE)
	out.usa[i,6] <- mean(temp$DV_A2[which(temp$RegimeA=="Democracy")],na.rm=TRUE)			
}

out.usa2 <- gather(as.data.frame(out.usa))
out.usa2$Outcome <- "Resolve in crises"
out.usa2$Outcome[grep("V4|V5|V6",out.usa2$key)] <- "Victory in war"
out.usa2$Regime <- "Dictatorship"
out.usa2$Regime[grep("V3|V6",out.usa2$key)] <- "Democracy"
out.usa2$Regime[grep("V2|V5",out.usa2$key)] <- "Transitional Democracy"
out.usa2$Regime <- as.factor(out.usa2$Regime)
out.usa2$Regime <- relevel(out.usa2$Regime, ref="Transitional Democracy")
out.usa2$Regime <- relevel(out.usa2$Regime, ref="Democracy")

out.plot <- ggplot(out.usa2, aes(y=value, fill=Regime)) + geom_density(alpha=0.5) + facet_wrap(~Outcome) + theme_bw() + ylab("Cell mean") + theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank()) + guides(fill=guide_legend(title="Regime Type"))
out.plot 